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Abstract 

A coupled finite element (FE) method and meshless 
local Petrov-Galerkin (MLPG) method for analyzing 
two-dimensional potential problems is presented in this 
paper. The analysis domain is subdivided into two 
regions, a finite element (FE) region and a meshless 
(MM) region. A single weighted residual form is 
written for the entire domain. Independent trial and test 
functions are assumed in the FE and MM regions. A 
transition region is created between the two regions. 
The transition region blends the trial and test functions 
of the FE and MM regions. The trial function blending 
is achieved using a technique similar to the ‘Coons 
patch’ method that is widely used in computer-aided 
geometric design. The test function blending is 
achieved by using either FE or MM test functions on 
the nodes in the transition element. 

The technique was evaluated by applying the coupled 
method to two potential problems governed by the 
Poisson equation. The coupled method passed all the 
patch test problems and gave accurate solutions for the 
problems studied. 

Keywords: Finite element method. Meshless Local 
Petrov-Galerkin (MLPG), Coupling method, Potential 
problem 

Introduction 

The finite element (FE) method [1, 2] is a widely used 
method, because of its accuracy, convenience, and 
flexibility. There are, however, deficiencies, such as 
element locking, discontinuous derivatives of the 
primary variable at element boundaries, and the need 
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for remeshing because of heavily distorted elements in 
large deformation. Meshless Methods (MM), such as 
the diffusion element method [3], the Element-Free 
Galerkin Method (EFGM) [4-6], h-p clouds [7], 
Partition of Unity (PU) [8], the Reproducing Kernel 
Particle Method (RKPM) [9], the Local Boundary 
Integral Equation (LBIE) method [10], and the 
Meshless Local Petrov-Galerkin (MLPG) method [11] 
are developed to avoid these deficiencies. Of these 
methods, the MLPG method is a very promising 
numerical method and has been successfully used on 
both potential and elasticity problems to obtain very 
accurate results [12-14]. The MLPG method is based on 
a local weak form and a moving least squares (MLS) 
approximation. This method is a “true" meshless 
method and does not need any “element" or “mesh” but 
uses a distributed set of nodes for both field 
interpolation and background integration. Although 
the MLPG method is attractive, the method is more 
computationally expensive than the FE method, and it 
is not as mature and comprehensive as the FE method. 

There is a great interest in combining two different 
numerical methods, because such coupling would 
exploit the potential of each method while avoiding 
their deficiencies. Several techniques have been 

proposed in the literature [15-22] to couple different 
numerical methods, such as combining the FEM with 
the Boundary Element Method (BEM) [15, 16], the 
FEM with the EFGM [17, 20-22], the EFGM with the 
BEM [18], the FEM with the MLPG method [19], etc. 
In coupling methods, the analysis domain is divided 
into two regions and each method is used in a region 
where it is most efficient. At the interface between the 
regions, the continuity of the primary variable and the 
balance of the secondary variables need to be 
maintained. One technique that is commonly used in 
the coupling methods [15-19] defines a transition 
region between the two regions. In the transition region, 
a ramp transition function is chosen to combine the trial 
functions (shape functions) of the two numerical 
methods. This technique is straightforward and widely 
used in coupling different numerical methods. Using 
this technique, the continuity condition of the primary 
variable is satisfied. However, the reciprocity' condition 
of the secondary variables is not guaranteed. Other 
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techniques, such as Lagrange multipliers method [20], a 
hierarchical mixed approximation method [21], and a 
method that comprises both FE and EFGM shape 
functions in the interface region [22] have also been 
presented recently. 

In this paper, a formulation that couples the FE and 
MLPG methods is presented for two-dimensional (2D) 
potential problems based on a single variational form. 
In this method, the analysis domain is divided into two 
regions (FEM region and MLPG region). Several 
computational issues related to the coupling method are 
discussed. Two numerical examples for the potential 
problems governed by Poisson equation are presented 
to evaluate the accuracy and efficiency of the proposed 
technique. 

Variational Coupling Procedure 

In this section the variational coupling procedure is 
described for coupling two disjoint regions (sub- 
domains) that are governed by a potential function. 

Consider a Poisson’s equation for a 2D problem 
governing the potential u (primary variable) in a two- 
dimensional (2D) domain Q bounded by contour f (see 
Figure 1) 

V 2 u = p (1) 

The right hand side p in Eq. (1) is, in general, a function 
of y and v. In the general case, the boundary T can have 
mixed boundary conditions. On one part of the 
boundary, T u , the potential u is prescribed, and on the 
remaining part T q , the secondary variables, flux, 
( q - du / dn ) are prescribed, /. e. 

u = u on T u and q = q on T q (2) 

where T = T u U T q . For example, in a steady-state heat 

transfer problem, which is a typical Poisson’s problem, 
u is a temperature function, q is the heat flux function 
and p is internal heat generation. 

The solution for Eq. (1) is sought in a weighted residual 
manner as 


J(|4 + |4-p)-v-<ft-c/y =0 (3) 

o dx ty 

where v = v(y,>>) is an arbitrary weight function. 

Using the divergence theorem, one can rewrite Eq. (3) 
as 


da dv du dv 
^ dx dx dy dy 

- J p • v • dQ = 0 


du 


C/OU cv OU < 7 V\ _ rO U 

“fe^ + ^^ ) ^ Q + r fe^ + ^ M ' ) - V£/r 


(4) 


where n x and n v are the direction cosines of the outward 
normal n of the boundary in the y and y directions, 
respectively (see Figure 1). 


In the coupling method, the domain Q is divided into an 
FE region (Q F ) and an MM region (Q M ) as shown in 
Figure 2. The subscripts F and M are used to describe 
the variables in the FE and MM regions, respectively. 
The interface between the two regions is denoted as /]. 

In the coupling method, the weighted residual equation 
(Eq. (3)) and subsequently the weak form (Eq. (4)) are 
written for the FE and MM regions separately as 


r { du dv du dv r du du 

- (—•- + - )'dQ + ( n x + » )-v dT 

dx ox dy dy ^ dr dv 


dy 


+ J(| ~ n i-; + ~»,.y)-v-dr - J p v dn=o 


dy 


n,. 


(5) 


K: = 


f ,du dv du dv r ou Ou 


c , du OU v r 

+ J C— - dT - J p v dn= 0 . 


du 




( 6 ) 

Notice that the outward normals of interface boundary 
from FE and MM domains are in opposite directions 
(see Figure 2), i.e. 


"i- = -n M (7) 

Therefore, the terms along the interface boundary' 
should meet the requirement 


r.du ou 


r - dx 

1 1 




dx 


du 

du 

dy 


( 8 ) 


n Mv )vdr = 0 


The balance of the secondary variables represented by 
Eq. (8) along the interface needs to be satisfied. This 
issue will be discussed in detail later in this paper. 
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r .8u dv du d\\ r f ou ou 

- [ ( + ) dQ. + (— -n x + — *n v )*v* 

t dx dx dy dy - dx dv 


r i' 

,du dv du 


dY 


p v dn - ( +— *— 

{ J rJ dx dx d\' dy> 

ii/,' i 2 A/ 

+ f (— n +— n ) v dr- f pvdn =0 

J rV dv ■ J 

L, 0X °> fiw 


(9) 


The next step in the coupling method is to choose the 
trial and test functions. The approximations for u in the 
variational form (Eq. (9)) are called the trial functions 
while the weight functions v in the variational form are 
called the test functions. The trial and test functions of 
the FE and MM methods along the interface (T|) are not 
the same. The coupling between the FE and MM 
regions can be achieved by defining a transition region 

Q t (see in Figure 3.) The region Q F ' is split into two 
parts, Q f and Q T with an interface T t . The region Q T is 
now between the FE region Q F and the mesh less region 


Trial and test functions 


The trial and test functions for the FE and MM regions 
are presented first and then the functions in the 
transition region are discussed. The trial and test 
functions are chosen such that 


u = Uf, + u M and v = v,. + v M (1 0) 


The trial functions of a point x in the FE region are 


u r 



0 


xe 

*€ Q m 


( 11 ) 


where u [ p are the interpolation functions used in the 

traditional FE method, and N f . is the number of the 
elements in the FE region. Similarly, the trial functions 
in the MM region are 


u 


M 


0 x e Qj, 

J a 'a/ 


( 12 ) 


where N M is the number of nodes in the MM region and 
(w A , ) ; are the moving least squares (MLS) 

approximations at each of the nodes in the MM region. 
The test functions are chosen as 





xe 

e=\ 


0 

*6 


(13) 


where are chosen to be the variations on u) c) (i.e. 
Su ¥ ) as is customary in the traditional FE method and 


0 


xe Q F 


( v a/ ) i x e 

l/=> 


(14) 


where (v A/ ), are chosen to be the Petrov-Galerkin weight 
function on a compact support sub-domain of each node 
[13]. The weight function of node / equals unity at x t , 

and monotonically decreases as || jc — || increases and 

is equal to zero at the boundary and outside of the sub- 
domain. 


Trial functions in the transition region 

The trial functions in the transition region need special 
attention. Different numerical methods (MM, BEM, 
and different types of elements of FEM) use distinct 
trial functions. Directly coupling two numerical 
methods by only matching the interface nodal values 
may not satisfy the compatibility condition along the 
interface. In other words, the function values produced 
in different regions along the interface may be different 
although a few nodal values are the same. Therefore, a 
transition element between the two regions that 
smoothly transforms the values of the trial functions 
from one region to another is preferred. 

A blending technique is developed to generate the 
interpolation function in a 4-sided region (quadrilateral 
element) in terms of its edge (boundary) functions (see 
Figure 4). The edges of the element can be curved and 
the number of the nodes on each of the edges need not 
be the same. This technique is similar to a technique, 
‘Coons patch’ [23], that is used in computer-aided 
geometric design and FEM mesh generation. The 
Coons patch technique is used to construct a smooth 
surface using the bounded curves. By using the Coons 
patch technique in an analogous fashion, the various 
trial functions on the edges of the transition elements 
are blended smoothly. This technique is also similar to 
the method used to generate serendipity elements. 
However, instead of using the boundary nodal values to 
generate the serendipity element, the boundary 
functions are used in this blending technique to 
construct the transition element. These boundary 
functions may be explicit functions of the boundary 
nodal values, such as FEM, or implicit functions of the 
nodal values, such as MM. 
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A bilinear blending is used in this study to construct the 
trial function of the element. Although the blending is 
bilinear , the trial functions can be an arbitrary ? function 
that depends on the edge functions of the element [23]. 
To be compatible with the FE method, the trial 
functions of the transition element are developed in a 
square element (parent element) in the 7* domain, 
where -I <(£, ;;)< 1 . The trial function can be written 
as 

“=j(i-£)-rj(7)+j(i+£)-r 4 (7) 

+i(i-i7)-r l (£)+l(i+7).r J (£) 

. “ i (,5) 


where r t (/ = 1,4) are the boundary functions of the 
quadrilateral element and Uj(j = n, p , q, and m) are the 
values of the primary variable at the comer nodes. 

The derivatives of the primary variable with respect to a 
and v can be obtained by calculating the derivatives 
with respect to £ and /;, then forming the Jacobian 
matrix of the element and following standard FE 
procedures. 

Typical transition elements between the FE region (Q F ) 
and MM region (Q M ) are shown in Figure 5. The open 
circles and solid circles represent the nodes in the FE 
and MM regions, respectively. The shaded elements are 
the transition elements. In these transition elements, the 
solid lines represent the boundaries that are connected 
to finite elements. On these boundaries, the edge 
functions can be generated by FE trial functions. The 
dashed lines in Figure 5 are the boundaries that connect 
to the MM region, and the edge functions are generated 
by the MM trial functions. The next step is to form the 
test function for the transition elements. 

Test functions in the transition region 

Although the test function can be arbitrarily selected in 
the weighted residual method, to satisfy the weak form 
requirement of the reciprocity conditions of the 
secondary variables with constant derivatives (see Eq. 
(8)), the test functions v must (a) be the same function 
and (b) have the same integration length along the 
interface of both of the regions. Therefore, the test 
function of a node in the transition element depends on 
the adjacent regions of that node. If the node is adjacent 
only to the FE region, and not to the MM region, the FE 
test function is used and the integration length of the 
function is identical to the side of the element (for 


example, the side lengths are a and b in type C of 
Figure 5 for node g). For a node that connects to the 
MM region, the MM test function with a support 
domain of /» is used (for example nodes q and k of type 
B in Figure 5). 

One additional concern about the test function of the 
transition region is the size of the support domain of a 
node using an MM test function. The size of the support 
domain of an MM node can be defined by a process 
similar to that used in references [11-14]. Although 
there is no restriction on the size of the support domain 
by the weighted residual method, the radius (/ 0 ) of the 
support domain of a node is defined to be less than the 
size of the element to which that node belongs (see for 
example, nodes q and k of type B in Figure 5). This 
limitation on the size of the support domain eliminates 
complicated numerical integration that may be needed 
to integrate accurately the weak form of the transition 
element. 

In Figure 5, three types of possible transition elements, 
A, B, and C, are identified. In the type B element, the 
nodes k and q are on the interface Tj. Therefore, the 
MM trial function is applied on edge kq , and the MM 
test functions should be applied on both nodes k and q. 
Note the sizes of the support domain (/„) of the test 
function for nodes k and q need not be the same. Since 
node q is shared by elements of type A and B, in order 
to satisfy the reciprocity conditions of the secondary 
variables on the common edge pq , the MM test function 
should be applied on node q of the element of type A, 
although no MM trial function is applied on any edge in 
element of type A. 

System equations 

The trial and test functions outlined in Eqs. (10-15) are 
substituted into the weak form in Eq. (9) and the 
integrations are performed. This leads to a system of 
equations similar to the FE and MLPG methods 

[KHu} = m (16) 

where [K] contains the [K F ], [K M ] and [K T ] matrices 
corresponding to the FE, MM and transition regions, 
respectively. Similarly, {f} contains contributions from 
the three regions. The solution of the system of 
equations leads to the potentials {u} at all the nodes in 
the model. Note that instead of working with fictitious 
nodal potentials {£} in the MM method the nodal 
potentials {w} are used by making use of the procedures 
contained in reference [14]. 

Boundary conditions 

There are no special considerations necessary to apply 
essential boundary conditions on the model. However, 
the values of the natural boundary condition (‘force’) 
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depend on the type of nodes (FE or MM) on the 
transition element. Based on Eq.(4), the 'force’ value 
on a node is calculated as 

/,= k~ ■",+—">) v rfr - \p \ (> 7 ) 

r, dx & a, 

where v f is the test function of the nodes j , F f is the 
boundary of the test function where the secondary 
variables, fluxes, are prescribed, and Qj is the domain 
of the test function. 

The formulation presented above for the coupled 
method is evaluated by performing patch test problems 
and by application to problems for which exact 
solutions are available. 

Numerical Examples 

To evaluate the algorithm of coupling FE-MM for 
potential problems, two example problems governed by 
Poisson’s equation, as shown in Figure 6, are 
considered. For the numerical integration in the MM 
domain, 8 and 12 Gauss points are used in the radial 
and tangential directions, respectively, of the circular 
support domain. A 3x3 quadrature rule is used for the 
quadratic elements in the FE region. The Gauss 
quadrature of a node in the transition element is 
dependent on the test function of the node. For the node 
using the MM test function, the same Gauss quadrature 
is used as in the MM region. A 4x4 quadrature rule is 
used for the node using the FE test functions. 

The potential and its partial derivatives obtained using 
the coupling method are compared to exact solutions 
wherever available. The error norm (jl^A/lb) is used to 
evaluate the effectiveness of various parameters. This 
norm is defined as 


where M is the total number of randomly distributed 
internal points in the domain at which the numerical 
solution is evaluated and compared to the exact 
solution. Note that these internal points are independent 
points that are randomly generated and are not 
associated with the nodes used in the models. In the 
MM region, a value of M = 50 is used, and in the FE 
region, six randomly distributed points are generated 
for each element. The error norm is calculated for the 
entire domain. If the exact solution of the problem is 
not available, the results of the coupling method are 
compared to the full FE solutions and the full MM 
solutions over the analysis domain. 

Two types of interfaces are used to model the examples. 
The first type of interface cuts across the domain and 
intersects the domain boundary as shown in Figure 2(a). 


The second type of interface forms a closed loop and 
does not intersect the domain boundary (see Figure 
2(b)). Two types of models are used to evaluate this 
type of interface. The first model uses an FE region that 
is embedded in an MM region, and the second model 
uses an MM region that is embedded in an FE region. 

In the examples (see Figure 6), 8-node quadratic 
isoparametric elements and quadratic basis functions in 
the MM region are used. When the coupled method 
reproduces the exact solution, then the method “passes" 
the patch tests. 

Example 1 

A potential problem governed by the Poisson’s 
equation, Vw 2 = 2 -(c, + c 2 ) is considered in a 
rectangular domain (Figure 6(a)). The exact solution for 
the problem is w = c,jc 2 + c 2 y 2 +c 3 xy, where c , are 
arbitrary constants. All four models in Figure 7 are used 
to evaluate the coupling method using quadratic 
quadrilateral elements and a quadratic basis function in 
the FE and MM regions, respectively. A type-I interface 
is used in the first model, and a type-1 1 interface is used 
in the models 2, 3 and 4. In these patch tests, the 
potential u is prescribed on the boundaries of the entire 
domain corresponding to the exact values. As expected, 
the exact solutions for both potential and fluxes in the 
entire domain are reproduced to machine accuracy. 
Similar results are obtained using the mixed boundary 
conditions. 


Example 2 

The second example considered is that of the torsion of 
a shaft of elliptical cross-section as shown in Figure 
6(b). The governing equation is a Poisson’s equation in 
terms of Prandtl stress function u with the boundary 
condition « = 0 on the boundary of the cross-section. 
The shear stresses in the shaft are determined as the 
derivatives of the stress function. For this problem, the 
right-hand side function of Eq. (1) is defined as 


p(x,y) = -2G0 = - 


2-A/V+6 2 ) 


(19) 


where G is the modulus of rigidity, 6 is the twist angle, 
K1 is the torsion moment applied on the shaft, and a and 
b are the semi-major and semi-minor axes of the elliptic 
cross-section, respectively. The exact solution of the 
problem [24] is 


u(x,y) = - 


M 


r— y2 n 

7z-a'h a 2 b 2 


( 20 ) 


Two coupling models are used for the problem, a full 
model of the elliptic shaft and a model of one-quarter of 
the shaft, as shown in Figure 8. As the exact solution is 
a quadratic, quadratic elements and quadratic basis 
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functions are used. In the full model, 16 quadratic 
elements (8-node quadrilateral) and 102 MM nodes 
using a quadratic basis function are used. The essential 
boundary condition of u - 0 is applied on the outer 
surface of the elliptic cross-section. In the model of 
one-quarter of the elliptic shaft, 4 quadratic elements 
and 38 nodes are used. The essential boundary 
condition (w = 0) is applied on the outer surface and 
the natural boundary conditions (</-0) are applied on 
the symmetric centerlines (y = 0 and x = 0). 

The error norm of the results for the potential and its 
first derivatives are computed using the exact solution 
at 100 randomly distributed points in the MM region 
and 96 random points in the FE region for the whole 
model. For the model of one-quarter of the elliptic 
shaft, 50 and 24 points are used in the MM and FE 
regions, respectively. The exact solutions are 
reproduced to machine accuracy for both models. 

Concluding Remarks 

A technique is presented for coupling the traditional 
finite element (FE) method and the meshless method 
(MM) for analyzing two-dimensional potential 
problems that are governed by Poisson’s equation. The 
coupling method allows the domain to be analyzed by 
the FE method in one region and the MM in another 
region. The regions are chosen depending on where 
each of methods is most effective. 

A single variational weighted residual formulation is 
used on the entire domain to develop the coupling 
method. A transition region is defined between the FE 
region and the MM region. In the transition region, a 
blending technique is used to construct the trial function 
based on the trial function in the adjacent region. The 
MM and FE test functions are used in the transition 
region to satisfy the compatibility condition of the 
secondary variables. 

Potential patch test problems involving Poisson 
equation are used to evaluate the effectiveness of the 
coupling technique when the meshless local Petrov- 
Galerkin method is used in the MM region. The 
coupling method passed all the patch test problems and 
yielded very accurate solutions for the problems 
studied. 
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(a) Type-I interface (b) Type-II interface 

Figure 2 Finite element and meshless method regions 



Figure 3 Finite element (FE) region, transition (TR) region and 

meshless method (MM) region 
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Region 2 



(a) Transition element 
in ‘x, y’ domain 



(b) Parent element 
in^ri’ 


Figure 4 Construction of transition element 



Figure 5 Typical transition elements 
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(a) Whole elliptic shaft model 



(b) Model of one-quarter of elliptic shaft 


Figure 8 Coupling models used in the elliptic shaft in torsion 
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Abstract 

A coupled finite element (FE) method and meshless 
local Petrov-Galerkin (MLPG) method for analyzing 
two-dimensional potential problems is presented in this 
paper. The analysis domain is subdivided into two 
regions, a finite element (FE) region and a meshless 
(MM) region. A single weighted residual form is 
written for the entire domain. Independent trial and test 
functions are assumed in the FE and MM regions. A 
transition region is created between the two regions. 
The transition region blends the trial and test functions 
of the FE and MM regions. The trial function blending 
is achieved using a technique similar to the ‘Coons 
patch’ method that is widely used in computer-aided 
geometric design. The test function blending is 
achieved by using either FE or MM test functions on 
the nodes in the transition element. 

The technique was evaluated by applying the coupled 
method to two potential problems governed by the 
Poisson equation. The coupled method passed all the 
patch test problems and gave accurate solutions for the 
problems studied. 

Keywords: Finite element method. Meshless Local 
Petrov-Galerkin (MLPG), Coupling method. Potential 
problem 

Introduction 

The finite element (FE) method [1,2] is a widely used 
method, because of its accuracy, convenience, and 
flexibility. There are, however, deficiencies, such as 
element locking, discontinuous derivatives of the 
primary variable at element boundaries, and the need 
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for remeshing because of heavily distorted elements in 
large deformation. Meshless Methods (MM), such as 
the diffusion element method [3], the Element-Free 
Galerkin Method (EFGM) [4-6], h-p clouds [7], 
Partition of Unity (PU) [8], the Reproducing Kernel 
Particle Method (RKPM) [9], the Local Boundary 
Integral Equation (LBIE) method [10], and the 
Meshless Local Petrov-Galerkin (MLPG) method [11] 
are developed to avoid these deficiencies. Of these 
methods, the MLPG method is a very promising 
numerical method and has been successfully used on 
both potential and elasticity problems to obtain very 
accurate results [12-14]. The MLPG method is based on 
a local weak form and a moving least squares (MLS) 
approximation. This method is a “true" meshless 
method and does not need any “element" or “mesh" but 
uses a distributed set of nodes for both field 
interpolation and background integration. Although 
the MLPG method is attractive, the method is more 
computationally expensive than the FE method, and it 
is not as mature and comprehensive as the FE method. 

There is a great interest in combining two different 
numerical methods, because such coupling would 
exploit the potential of each method while avoiding 
their deficiencies. Several techniques have been 
proposed in the literature [15-22] to couple different 
numerical methods, such as combining the FEM with 
the Boundary Element Method (BEM) [15, 16], the 
FEM with the EFGM [17, 20-22], the EFGM with the 
BEM [18], the FEM with the MLPG method [19], etc. 
In coupling methods, the analysis domain is divided 
into two regions and each method is used in a region 
where it is most efficient. At the interface between the 
regions, the continuity of the primary variable and the 
balance of the secondary variables need to be 
maintained. One technique that is commonly used in 
the coupling methods [15-19] defines a transition 
region between the two regions. In the transition region, 
a ramp transition function is chosen to combine the trial 
functions (shape functions) of the two numerical 
methods. This technique is straightforward and widely 
used in coupling different numerical methods. Using 
this technique, the continuity condition of the primary' 
variable is satisfied. However, the reciprocity condition 
of the secondary variables is not guaranteed. Other 
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techniques, such as Lagrange multipliers method [20], a 
hierarchical mixed approximation method [21], and a 
method that comprises both FE and EFGM shape 
functions in the interface region [22] have also been 
presented recently. 

In this paper, a formulation that couples the FE and 
MLPG methods is presented for two-dimensional (2D) 
potential problems based on a single variational form. 
In this method, the analysis domain is divided into two 
regions (FEM region and MLPG region). Several 
computational issues related to the coupling method are 
discussed. Two numerical examples for the potential 
problems governed by Poisson equation are presented 
to evaluate the accuracy and efficiency of the proposed 
technique. 

Variational Coupling Procedure 

In this section the variational coupling procedure is 
described for coupling two disjoint regions (sub- 
domains) that are governed by a potential function. 

Consider a Poisson’s equation for a 2D problem 
governing the potential u (primary variable) in a two- 
dimensional (2D) domain Q bounded by contour T (see 
Figure 1) 


V 2 u = p (1) 

The right hand side p in Eq. (1) is, in general, a function 
ofx an dy. In the general case, the boundary T can have 
mixed boundary conditions. On one part of the 
boundary, T u , the potential u is prescribed, and on the 
remaining part F q , the secondary variables, flux, 
( q ~ dui dn ) are prescribed, i. e. 

u = u on T u and q = q on T q (2) 

where T = T u U T q . For example, in a steady-state heat 

transfer problem, which is a typical Poisson’s problem, 
w is a temperature function, q is the heat flux function 
and p is internal heat generation. 

The solution for Eq. (1) is sought in a weighted residual 
manner as 


c,d 2 u d 2 u , 

+ — 7-p)-v-dx-dy =0 ( 3 ) 

n dx $> r 

where v = v(x,y) is an arbitrary weight function. 

Using the divergence theorem, one can rewrite Eq. (3) 
as 


c du dv du dv rdu 

- (— *— + )</q + ( — 

q dx dx dy dy f dx 

- J p-v- dQ = 0 


du 

+ n )-v dY 

dy ' 


(4) 

where n x and n v are the direction cosines of the outward 
normal n of the boundary in the x and y directions, 
respectively (see Figure 1). 


In the coupling method, the domain Q is divided into an 

FE region (Q F ) and an MM region (Q M ) as shown in 
Figure 2. The subscripts F and M are used to describe 
the variables in the FE and MM regions, respectively. 
The interface between the two regions is denoted as /}. 

In the coupling method, the weighted residual equation 
(Eq. (3)) and subsequently the weak form (Eq. (4)) are 
written for the FE and MM regions separately as 


// = 


r du dv du dv r m du 

J + J^- n > + ^ "^ v dT 

du 

dy 


n ',dxdx dy dy' \]S ^ ' fy 

r ,du du x r 

+ r J ( &'" ,<+ ^' w '-' ) ' v ' i/r ~J p- v dQ =° 


n 


r' 


(5) 


h; 


r du dv du ov rdu du 

~ J + 1 (tt ”. +— ' rt v)' v dr 

n J dx fa dy dy r *dx dy > 

M AY 

r ,du du „ „ r 

+ J + ^T'"m.) v dT - J p v dCl =0 . 

( 6 ) 


Notice that the outward normals of interface boundary 
from FE and MM domains are in opposite directions 
(see Figure 2), i.e. 


= -n„ (7) 

Therefore, the terms along the interface boundary 
should meet the requirement 


c ,du du x 
r ,du du . 

v ( & 


( 8 ) 


The balance of the secondary variables represented by 
Eq. (8) along the interface needs to be satisfied. This 
issue will be discussed in detail later in this paper. 


2 

American Institute of Aeronautics and Astronautics 


AIAA-2002-1659 


Combing Eqs. (5) and (6) and using Eq. (8), a single 
variational equation that is valid for the entire domain 
Li can be written as 


v,. 



jce fi F 


(13) 


_ ( ^L.^^yda + l 

dx dx dy dy jy, dx dy 


dV 




r r y du dv du dx\ 

j p-v-dQ - j 

J rT d X OX ay d\’ 


■b’ 


+ f n v )vdr- [ p-vdQ =0 

r > dy n, 


(9) 


The next step in the coupling method is to choose the 
trial and test functions. The approximations for u in the 
variational form (Eq. (9)) are called the trial functions 
while the weight functions v in the variational form are 
called the test functions. The trial and test functions of 
the FE and MM methods along the interface (Tj) are not 
the same. The coupling between the FE and MM 
regions can be achieved by defining a transition region 

Q t (see in Figure 3.) The region Q F ' is split into two 
parts, Qf and Q T with an interface T t . The region Q T is 
now between the FE region Q F and the meshless region 
Qm. 


Trial and test functions 


The trial and test functions for the FE and MM regions 
are presented first and then the functions in the 
transition region are discussed. The trial and test 
functions are chosen such that 


u = u } + u M and v = iy +v M (1 0) 

The trial functions of a point x in the FE region are 


x e n F 

( >= l 

0 X€ Q., 


(ID 


where mJ. 1 ’ are the interpolation functions used in the 
traditional FE method, and N,. is the number of the 
elements in the FE region. Similarly, the trial functions 
in the MM region are 


0 

£(»«), 

7=1 


x e Q, 

x e n M 


( 12 ) 


where N M is the number of nodes in the MM region and 
(« A/ ) are the moving least squares (MLS) 

approximations at each of the nodes in the MM region. 
The test functions are chosen as 


where v^are chosen to be the variations on (/.*?. 
du ? ) as is customary in the traditional FE method and 


0 

2 > m ), 

7=1 


xe Q F 
X ^ 


(14) 


where (v A/ ) y are chosen to be the Petrov-Galerkin weight 
function on a compact support sub-domain of each node 
[13]. The weight function of node / equals unity at x t , 

and monotonically decreases as ||jc-^|| increases and 
is equal to zero at the boundary and outside of the sub- 
domain. 

Trial functions in the transition region 

The trial functions in the transition region need special 
attention. Different numerical methods (MM, BEM, 
and different types of elements of FEM) use distinct 
trial functions. Directly coupling two numerical 
methods by only matching the interface nodal values 
may not satisfy the compatibility condition along the 
interface. In other words, the function values produced 
in different regions along the interface may be different 
although a few nodal values are the same. Therefore, a 
transition element between the two regions that 
smoothly transforms the values of the trial functions 
from one region to another is preferred. 

A blending technique is developed to generate the 
interpolation function in a 4-sided region (quadrilateral 
element) in terms of its edge (boundary) functions (see 
Figure 4). The edges of the element can be curved and 
the number of the nodes on each of the edges need not 
be the same. This technique is similar to a technique, 
‘Coons patch’ [23], that is used in computer-aided 
geometric design and FEM mesh generation. The 
Coons patch technique is used to construct a smooth 
surface using the bounded curves. By using the Coons 
patch technique in an analogous fashion, the various 
trial functions on the edges of the transition elements 
are blended smoothly. This technique is also similar to 
the method used to generate serendipity elements. 
However, instead of using the boundary nodal values to 
generate the serendipity element, the boundary 
functions are used in this blending technique to 
construct the transition element. These boundary 
functions may be explicit functions of the boundary 
nodal values, such as FEM, or implicit functions of the 
nodal values, such as MM. 
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A bilinear blending is used in this study to construct the 
trial function of the element. Although the blending is 
bilinear , the trial functions can be an arbitrary’ function 
that depends on the edge functions of the element [23], 
To be compatible with the FE method, the trial 
functions of the transition element are developed in a 
square element (parent element) in the % rf domain, 
where -I <(£, rf)< I . The trial function can be written 
as 

« = ^(i-^)-r 3 (/7)+^(i+4 : )r 4 (;7) 

+^(i->7)-r l (£)+^(i+t7).r J (£) 

. , (,5) 

~(1 + £Xl + >?)•«, "(I -£XI+7)-M», 


where r, (/ = 1,4) are the boundary functions of the 
quadrilateral element and u, (j = n, p , q, and m) are the 
values of the primary variable at the comer nodes. 

The derivatives of the primary variable with respect to x 
and y can be obtained by calculating the derivatives 
with respect to £ and //, then forming the Jacobian 
matrix of the element and following standard FE 
procedures. 

Typical transition elements between the FE region (Q F ) 
and MM region (Q M ) are shown in Figure 5. The open 
circles and solid circles represent the nodes in the FE 
and MM regions, respectively. The shaded elements are 
the transition elements. In these transition elements, the 
solid lines represent the boundaries that are connected 
to finite elements. On these boundaries, the edge 
functions can be generated by FE trial functions. The 
dashed lines in Figure 5 are the boundaries that connect 
to the MM region, and the edge functions are generated 
by the MM trial functions. The next step is to form the 
test function for the transition elements. 

Test functions in the transition region 

Although the test function can be arbitrarily selected in 
the weighted residual method, to satisfy the weak form 
requirement of the reciprocity conditions of the 
secondary variables with constant derivatives (see Eq. 
(8)), the test functions v must (a) be the same function 
and (b) have the same integration length along the 
interface of both of the regions. Therefore, the test 
function of a node in the transition element depends on 
the adjacent regions of that node. If the node is adjacent 
only to the FE region, and not to the MM region, the FE 
test function is used and the integration length of the 
function is identical to the side of the element (for 


example, the side lengths are a and h in type C of 
Figure 5 for node g). For a node that connects to the 
MM region, the MM test function with a support 
domain of / 0 is used (for example nodes q and k of type 
B in Figure 5). 

One additional concern about the test function of the 
transition region is the size of the support domain of a 
node using an MM test function. The size of the support 
domain of an MM node can be defined by a process 
similar to that used in references [11-14]. Although 
there is no restriction on the size of the support domain 
by the weighted residual method, the radius (/ 0 ) of the 
support domain of a node is defined to be less than the 
size of the element to which that node belongs (see for 
example, nodes q and k of type B in Figure 5). This 
limitation on the size of the support domain eliminates 
complicated numerical integration that may be needed 
to integrate accurately the weak form of the transition 
element. 

In Figure 5, three types of possible transition elements, 
A, B, and C, are identified. In the type B element, the 
nodes k and q are on the interface rj. Therefore, the 
MM trial function is applied on edge kq , and the MM 
test functions should be applied on both nodes k and q. 
Note the sizes of the support domain (/,,) of the test 
function for nodes k and q need not be the same. Since 
node q is shared by elements of type A and B, in order 
to satisfy the reciprocity conditions of the secondary 
variables on the common edge pq , the MM test function 
should be applied on node q of the element of type A, 
although no MM trial function is applied on any edge in 
element of type A. 

System equations 

The trial and test functions outlined in Eqs. (10-15) are 
substituted into the weak form in Eq. (9) and the 
integrations are performed. This leads to a system of 
equations similar to the FE and MLPG methods 

[K]-{u} = {f} (16) 

where [K] contains the [K F ], [K M ] and [K T ] matrices 
corresponding to the FE, MM and transition regions, 
respectively. Similarly, {f} contains contributions from 
the three regions. The solution of the system of 
equations leads to the potentials {u} at all the nodes in 
the model. Note that instead of working with fictitious 
nodal potentials {£} in the MM method the nodal 
potentials {u} are used by making use of the procedures 
contained in reference [14]. 

Boundary conditions 

There are no special considerations necessary to apply 
essential boundary conditions on the model. However, 
the values of the natural boundary condition (‘force') 
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depend on the type of nodes (FE or MM) on the 
transition element. Based on Eq.(4), the ‘force' value 
on a node is calculated as 

/,-/<! -Ip-*, ■<& O’) 

r J & 8y n 

where v, is the test function of the nodes j , r f is the 
boundary of the test function where the secondary 
variables, fluxes, are prescribed, and Q, is the domain 
of the test function. 

The formulation presented above for the coupled 
method is evaluated by performing patch test problems 
and by application to problems for which exact 
solutions are available. 

Numerical Examples 

To evaluate the algorithm of coupling FE-MM for 
potential problems, two example problems governed by 
Poisson’s equation, as shown in Figure 6, are 
considered. For the numerical integration in the MM 
domain, 8 and 12 Gauss points are used in the radial 
and tangential directions, respectively, of the circular 
support domain. A 3x3 quadrature rule is used for the 
quadratic elements in the FE region. The Gauss 
quadrature of a node in the transition element is 
dependent on the test function of the node. For the node 
using the MM test function, the same Gauss quadrature 
is used as in the MM region. A 4x4 quadrature rule is 
used for the node using the FE test functions. 

The potential and its partial derivatives obtained using 
the coupling method are compared to exact solutions 
wherever available. The error norm (||e M || 2 ) is used to 
evaluate the effectiveness of various parameters. This 
norm is defined as 


h I, - Ji>, Jj M < 1 8 > 

where M is the total number of randomly distributed 
internal points in the domain at which the numerical 
solution is evaluated and compared to the exact 
solution. Note that these internal points are independent 
points that are randomly generated and are not 
associated with the nodes used in the models. In the 
MM region, a value of M = 50 is used, and in the FE 
region, six randomly distributed points are generated 
for each element. The error norm is calculated for the 
entire domain. If the exact solution of the problem is 
not available, the results of the coupling method are 
compared to the full FE solutions and the full MM 
solutions over the analysis domain. 

Two types of interfaces are used to model the examples. 
The first type of interface cuts across the domain and 
intersects the domain boundary as shown in Figure 2(a). 


The second type of interface forms a closed loop and 
does not intersect the domain boundary (see Figure 
2(b)). Two types of models are used to evaluate this 
type of interface. The first model uses an FE region that 
is embedded in an MM region, and the second model 
uses an MM region that is embedded in an FE region. 

In the examples (see Figure 6), 8-node quadratic 
isoparametric elements and quadratic basis functions in 
the MM region are used. When the coupled method 
reproduces the exact solution, then the method “passes’ 
the patch tests. 

Example 1 

A potential problem governed by the Poisson s 
equation, Vw 2 =2*(c, + c 2 ) is considered in a 
rectangular domain (Figure 6(a)). The exact solution for 
the problem is u = c { x 2 + c 2 y 2 + c } xy , where c, are 
arbitrary constants. All four models in Figure 7 are used 
to evaluate the coupling method using quadratic 
quadrilateral elements and a quadratic basis function in 
the FE and MM regions, respectively. A type-1 interface 
is used in the first model, and a type-11 interface is used 
in the models 2, 3 and 4. In these patch tests, the 
potential u is prescribed on the boundaries of the entire 
domain corresponding to the exact values. As expected, 
the exact solutions for both potential and fluxes in the 
entire domain are reproduced to machine accuracy. 
Similar results are obtained using the mixed boundary 
conditions. 


Example 2 

The second example considered is that of the torsion of 
a shaft of elliptical cross-section as shown in Figure 
6(b). The governing equation is a Poisson's equation in 
terms of Prandtl stress function u with the boundary 
condition u = 0 on the boundary of the cross-section. 
The shear stresses in the shaft are determined as the 
derivatives of the stress function. For this problem, the 
right-hand side function of Eq. (1) is defined as 


p(x,y) = -2G6 = - 


2M -{cf+b 2 ) 
K'a * b } 


(19) 


where G is the modulus of rigidity, 6 is the twist angle, 
M is the torsion moment applied on the shaft, and a and 
b are the semi-major and semi-minor axes of the elliptic 
cross-section, respectively. The exact solution of the 
problem [24] is 


u(x,y) = 


M 


7t • a 


b 



•> 

r_ 

t r 


1) 


( 20 ) 


Two coupling models are used for the problem, a full 
model of the elliptic shaft and a model of one-quarter of 
the shaft, as shown in Figure 8. As the exact solution is 
a quadratic, quadratic elements and quadratic basis 
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functions are used. In the full model, 16 quadratic 
elements (8-node quadrilateral) and 102 MM nodes 
using a quadratic basis function are used. The essential 
boundary condition of u - 0 is applied on the outer 
surface of the elliptic cross-section. In the model of 
one-quarter of the elliptic shaft, 4 quadratic elements 
and 38 nodes are used. The essential boundary 
condition (// = 0) is applied on the outer surface and 
the natural boundary conditions (q = 0) are applied on 
the symmetric centerlines (y = 0 and x = 0). 

The error norm of the results for the potential and its 
first derivatives are computed using the exact solution 
at 100 randomly distributed points in the MM region 
and 96 random points in the FE region for the whole 
model. For the model of one-quarter of the elliptic 
shaft, 50 and 24 points are used in the MM and FE 
regions, respectively. The exact solutions are 
reproduced to machine accuracy for both models. 

Concluding Remarks 

A technique is presented for coupling the traditional 
finite element (FE) method and the meshless method 
(MM) for analyzing two-dimensional potential 
problems that are governed by Poisson’s equation. The 
coupling method allows the domain to be analyzed by 
the FE method in one region and the MM in another 
region. The regions are chosen depending on where 
each of methods is most effective. 

A single variational weighted residual formulation is 
used on the entire domain to develop the coupling 
method. A transition region is defined between the FE 
region and the MM region. In the transition region, a 
blending technique is used to construct the trial function 
based on the trial function in the adjacent region. The 
MM and FE test functions are used in the transition 
region to satisfy the compatibility condition of the 
secondary variables. 

Potential patch test problems involving Poisson 
equation are used to evaluate the effectiveness of the 
coupling technique when the meshless local Petrov- 
Galerkin method is used in the MM region. The 
coupling method passed all the patch test problems and 
yielded very accurate solutions for the problems 
studied. 
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Figure 1 Two-dimensional potential problem 
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(a) Type- 1 interface (b) Type-II interface 

Figure 2 Finite element and meshless method regions 



0 X 


Figure 3 Finite element (FE) region, transition (TR) region and 

meshless method (MM) region 
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Region 2 




(a) Transition element 
in ‘x, y’ domain 


(b) Parent element 
in T|’ 


Figure 4 Construction of transition element 



Figure 5 Typical transition elements 
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(a) Whole elliptic shaft model (b) Model of one-quarter of elliptic shaft 


Figure 8 Coupling models used in the elliptic shaft in torsion 
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